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Abstract 

In this note we apply a modified fractional Bessel differential equa- 
tion to the problem of describing corneal topography. We find the 
solution in terms of the power series. This solution has an interesting 
behavior at infinity which is a generalization of the classical results for 
modified Bessel function of order 0. Our model fits the real corneal 
geometry data with an error of order of a few per cent. 

1 Introduction 

Sight is the most crucial sense that we posses since it enables us to per- 
ceive the world very accurately. With the advance of medical technology 
treating various eye diseases becomes more adequate and successful. This 
would not be possible without proper mathematical models of biomechanics 
of eye and its constituents. One of the most important parts of the human 
eye is the cornea because it is responsible for about two-thirds of refractive 
power (for biological treatment see for ex. [lj). Mathematical description of 
corneal topography is very important from the point of view of ophthalmol- 
ogists because many seeing disorders originate in some distortions in corneal 
geometry. 

There are many types of corneal topography mathematical models. The 
most common and simple are based on conic sections (see for ex. [21 E]). 
Unfortunately, they are taken without much physical motivation. Some, 
more complex, models are based on finite element methods or shell theory 
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[4, 5j. Models that arc widely used to describe abberations in cornea or lens 
often use Zernike orthogonal polynomials j6j[7]. 

In [B] we have proposed yet another model of corneal topography. It was 
based on physical derivation and membrane equation. The first approxima- 
tion to that model is 

-Ah + ah = b, h\ m = 0, (1) 

where a and b are dimensionless, positive constants and Q is the domain 
on which cornea is situated. In this letter we generalize this model and use 
fractional derivatives instead of classical ones. As a result we obtain a modi- 
fied fractional Bessel differential equation which we solve and find interesting 
asymptotic behavior. The fractional calculus methods are increasingly more 
popular and very successful in a large number of physical application (see for 
ex. (9j dni E] ) • A comprehensive introductions to fractional calculus and its 
applications can be found for expample in [T2l [T3] . Up to authors' knowl- 
edge the fractional Bessel equation was investigated only in |14j . However, 
the form of that equation was different from the one analyzed by us. 



2 Model and its analysis 

Assume that the corneal surface is axisymmetric, that is h = h(r) where 
r G [0, 1] and Q is a unit circle. We also assume that h can be represented 
by convergent power series. Rewriting Q in polar coordinates we obtain 

--^-(rh') + ah = b, h'(0) = 0, h(l) = 0. (2) 
r dr 

The condition at r = grants us a smooth solution at the origin while 
h(l) = determines the rim of cornea. Now, multiply (12]) by r 2 , let y = — h 
and x = \J~ar to get 

x— (xy') - x 2 y = -x 2 } (3) 
dx a 

which is a nonhomogeneous modified Bessel equation of order 0. Its solution 
and application to corneal topography was presented in [8 J . Now, we want to 
generalize this model and use fractional derivative instead of the usual one. 
We propose 

x a D«{xy')-x 2 y= h -x 2 , y'(0) = 0, y(y/a) = 0, < a < 1, (4) 

as the model of corneal geometry. Operator Dq is the Riemann-Liouville 
fractional derivative defined by D£y := 4zDq a ^y, where 



(iV*v)(*):=r^ / {x-ty- l y{t)dt (5) 
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is the fractional integral operator (see |12j). It can be shown (see |T3| [16]) 
that in our setting since xy'(x)\ x= o = and y is analytic the Riemann-Lioville 
and Caputo fractional derivatives coincide. The choice of appropriate power 
of x in ^ is neccesary to maintain dimensional consistency and physical 
meaning. We will call Q the modified fractional Bessel equation of order 
in analogy with the classical case. 

Immediately we see that y p (x) = —b/a is the particular solution of Q. 
Because the general solution to Q has the form y = y^ + y p , we only have 
to find the solution yjj to the homogeneous equation 

x a D«(xy') = x 2 y, y'(0) = 0, < a < 1, (6) 

where the condition y'(0) = comes from the fact that y p is constant. We 
seek for a solution to ^ in the form of power series, that is 

oo 

y(x) = ^a n x n . (7) 

n=0 

Substituting it into @ and noting that (Dgt n )(x) = T(n+1) /T(n-a+l)x n - a 
we obtain recurrence formulas for coefficients a n 

r(n — a + 1) 

a n = -— -r a n _ 2 , n > 2. (8) 

n z l (n) 

The equations for do and a± are automatically fulfilled since = y'{0) = at 
and we can choose clq = 1. By the ratio test we see that the series is abso- 
lutely convergent. Thus, taking into account Q we obtain the solution to 
homogeneous fractional Bessel equation ^ which in accordance with classi- 
cal theory we will denote by Iq 



Vh{x) 

for some constant C and the convention that HiLi = 1- This modified 
fractional Bessel function is a generalization of the classical Bessel function 
since Jq = J . Additionally, we can see that 7° = exp(x 2 /2). Noting the 
other boundary condition y(y/a) = and returning to the original variables 
h and r we can write the solution to the boundary value problem Q 



This equation describes the shape of human cornea and we will see later that 
it gives very accurate fit with the real data. 
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The modified fractional Bessel function Iq defined in (|9]) has very inter- 
esting behavior as x — > oo. As we will see it is a generalization of the classical 
results from asymptotic theory. First, we prove a technical lemma. 

Lemma 1. Let F^x) := t^e~ xt dt, then has the following leading order 
behavior as x — > oo 

F M 0r)~x-^ +1 )r( M + l) (11) 



Proof. With a change of variable (s = xt) we can write the integral defining 



Ffj, as 



F^x) = x~^ +1) J s"e- s ds = x~^ +1) (r(fi + 1) - J s^e-'ds^j . (12) 
Now, integrating by parts we write the asymptotic expansion of F^ 

^.^r^ + D-.-gj^i^^, (is, 



as x — y oo. Since the second term in (13) is exponentially small the leading 



order behavior of F M is defined by te first term and that concludes the proof. 

□ 

The proof of following theorem is based on the Laplace method for asymp- 
totic integrals. The main point is that we are dealing with integro-differential 
equation which introduces some complications. 

Theorem 2. The solution to modified Bessel fractional differential equation 
(j^|) has the following asymptotic leading order behavior as x — > oo 

q(2-cp / 1 + a 2 \ 

Iq(x) ~ x !+« exp — - — . (14) 



Remark 1. We see that for a = 0, 1 asymptotic form (14) reduces to well 
known formulas for Iq(x) = exp(x 2 /2) and Iq(x) = Io(x) ~ x~ 1 ^ 2 exp(x) 
(see for ex. 115]). 

Proof. By using the identity for composition of fractional integral D^ a and 
Riemann-Lioville fractional derivative (|13j. formulas 2.113 and 2.135-136) 
and by the fact that xy'(x)\ x= o = we have 

(D a D«(xy\x)))(x)=xy\x) < a < 1. (15) 
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We transform Q into an Volterra integro-differential equation 
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y \x) 



Via) x 



(x - t) a -H 2 ' a y(t)dt. 



(16) 



Since we are looking for the behavior of ( 16 ) for large x we change the variable 



s = t/x to obtain constant limits of integration 

I {l-s) a ~ 1 s 2 - a y{sx)ds. 
Jo 



y [X) 



r (a) 



;i7) 



Now, guided by Remark [T] we seek the approximate solution of (16) in the 
form y(x) = f(x) exp (x x /\), where / and A are to be determined and / has 
algebraic growth. We have 



y [x) 



x 



T(a) 



'.1-8 



,a-l„2-a 



f(sx)e^ s ds. 



Notice that for large x the integrand in (18) is dominated by exponential 
term which has its maximum at s = 1. To move this maximum to the lower 
limit we substitute again t = 1 — s and obtain 



y (x) 



x 



Via) 



(19) 



As in the Laplace method for asymptotic integrals (see for ex. |17| ) for x — > 
oo the greatest contribution to the integral ( 16 ) comes from the neighborhood 
of the maximum of exponential term, that is when t is close to 0. Writing 
(1 — t) x = 1 — At + f(x(l — £)) = f(x) — xf'(x)t + ... and retrieving only 
first terms we can make an approximation valid for large x 



y [x) 



x 



r (a) 



-e * 



f(x) 



'\-ty- a t 



2 - a + a - l e~ xt dt 



+xf\x) (1 - tf- a t a e- xH d^j 



(20) 



In the first integral in (20) we further approximate (1 — t) 2 ~ a w 1 — (2 — a)t 
and in the second (1 — t) 2 ~ a w 1. This is justified since the whole mass of 
the integral is focused near t = 0. Using Lemma [T] we can write 



y'(x) « [f(x) (-F a -i(x A ) - (2 - a)F a (x x )) + xf(x)F a (x x )] 

~ [(a;- Aa+1 - a(2 - a)aT A ( Q+1 ) +1 ) f(x) - ax 2 ~ x{a+1) f '(x)] 



(21) 



where we used the formula T(a + 1) = aT(a). Noticing that y'(x) = 
f'{x) exp(x A /A) + x x ~ x exp(x A /A) we finally obtain a differential equation 
for / 

(1 + ax 2 - A(Q+1) ) /' = (x" Aq+1 - x*- 1 - a{2 - a)x~ x(a+1)+1 ) f. (22) 

We see that the choice of A = 2/(1 + a) is not only necessary for / to have 
algebraic growth but also greatly simplifies this equation into 

/' q(2-a) 1 

/ " 1 + a x' { ' 

which has the solution 

a(2 — a) 

f(x) = Cx~^+^. (24) 
This finishes the proof. □ 



3 Application to corneal topography 



In this section we apply previously analyzed mathematical model (10) to 
a dataset consisting of 123 x 123 measurement points of real cornea. We 
fit the solution h finding the least squares unknown parameters a, b and a. 
Cornea has its thickness thus we make two fits: one for exterior and one for 
interior surface. For exterior surface we find a = 0.580404, b = 1.19734 and 
a = 0.421345 with a mean absolute fitting error 0.014mm. Similarly, for 
interior surface we have a = 0.818763, b = 1.66664, a = 0.503431 and the 
mean absolute fitting error 0.03mm. We see that both parameters a and b 
are of order of unity what was expected (see [8]) and the order of derivative in 
(El) lies between 1 and 2. Contour plots of absolute fitting errors are depicted 
on Figure [T] 



4 Conclusion 

In this letter we have proposed a new mathematical model of corneal topog- 
raphy based on the modified fractional Bessel differential equation. We found 
its solution in a form of absolutely convergent power series. This solution is 
a generalization of the classical modified Bessel function of order 0. Also, we 
have found the behavior of Iq as x — > oo, which reduces to the known cases 
as a = 0, 1. When applying our model to the real corneal data we found 
accurate fit and parameter orders predicted by theory. 



6 



Absolute error for fitting (Exterior Surface) [mm] Absolute error for fitting (Interior Surface) [mm] 




Figure 1: Absolute fitting errors for exterior (left) and interior (right) sur- 
faces. 
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